Pocket Monte Carlo algorithm for classical doped dimer models 
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We study the correlations of classical hardcore dimer models doped with monomers by Monte 
Carlo simulation. We introduce an efficient cluster algorithm, which is applicable in any dimension, 
for different lattices and arbitrary doping. We use this algorithm for the dimer model on the square 
lattice, where a finite density of monomers destroys the critical confinement of the two-monomer 
problem. The monomers form a two-component plasma located in its high-temperature phase, 
with the Coulomb interaction screened at finite densities. On the triangular lattice, a single pair of 
monomers is not confined. The monomer correlations are extremely short-ranged and hardly change 
with doping. 
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I. INTRODUCTION 

Models of classical dimers on a lattice play an impor- 
tant role in polymer physics, but their statistical mechan- 
ics has proven to be of great value in a much wider range 
of settings. Prominently, in the early days of the resonat- 
ing valence bond theoriesu of high-temperature supercon- 
ductivity, Rokhsar and Kivelsona proposed that dimers 
could act as caricatures of singlet ('valence') bonds be- 
tween spins 1/2, with the hardcore condition on the dimer 
configurations translating the constraint that each spin 
participate in exactly one singlet bond with a neighboring 
spin. Resonances between different configurations gener- 
ate a quantum dynamics for the resulting Rokhsar Kivel- 
son quantum dimer models (RK-QDM). 

These RK-QDMs are also of interest as effective theo- 
ries, as they appear as limiting cases of d — 2 + 1 dimen- 
sional Ising sauge theories which are dual to frustrated 
Ising modelsfl Central in this context is the question of 
confinement - can two test charges be separated infinitely 
far at finite cost in free energy? In the dimer model, 
these test charges are monomers, representing spinons or 
holons in the RVB theory. 

In the context of RK quantum dimer models, classical 
dimer models play a crucial role because, for a partic- 
ular choice of parameters, the quantum wavefunction is 
an equal-amplitude superposition of all classical dimer 
configurations .□ For this reason, dimer-diagonal correla- 
tion functions can be obtained from the classical (equally 
weighted) average over all dimer configurations. At zero 
doping, such correlations can be obtained analytically us- 
ing Pfaffian methods introduced by KastjcleynQ and de- 
veloped further by a number of authors!] These meth- 
ods have been extended to obtain correlation, functions 
of a pair of monomers on the square lattice,Lrcj although 
it is not clear whether asymptotic monomer correlations 
can in general be obtained in closed form along those 
lines. Present analytical approaches do not allow an 
exact treatment of finite monomer densities; this is re- 
lated to the absence of an analytical solution, of the two- 
dimensional Ising model in a magnetic field.u 

To reach the RK point in the presence of doping, it 



is necessary not only to fine-tune two ratios of magni- 
tudes of kinetic and potential energies but also a hopping 
phaseJj This procedure is not innocuous, but classical 
correlations are also interesting because they are those of 
the quantum model at high temperature. There, hopping 
and resonance as well as any potential energies included 
in the dimer model are much smaller than the thermal 
energy. The classical correlations can remain non-trivial 
as a result of the projection onto the space of monomer- 
dimer coverings. One is. thus interested in the presence 
of (at least short-rangeEJ) correlated structures, and how 
they change with doping. 

In this paper, we present an efficient numerical al- 
gorithm for simulating classical monomer-dimer models. 
The algorithm is analyzed in some detail, and we mention 
possible uses for related problems. 

Our main results are the following. Firstly, a pair of 
monomers on the triangular lattice is deconfined, with a 
correlation length of less than one lattice constant. This 
is in keeping with the very short-range dimer correla- 
tions on that latticje.liil On the square lattice, the critical 
dimer correlations!] generate a Coulomb interaction, with 
monomers on the same (different) sublattice having equal 
(opposite) charges. At finite density, the monomers form 
a two-dimensional two-component plasma in its Debye- 
screened high-temperature phase. For both square and 
triangular lattice, the monomer correlations decay mono- 
tonically with increasing density, without any sign of ad- 
ditional new correlations on any length scale. 

We note that the^tudy of monomer-dimer models goes 
back a long time.Ej However, we have not found in the 
literature an evaluation of the correlations studied here. 



II. POCKET ALGORITHM FOR MONTE 
CARLO SIMULATION OF DIMER MODELS 

Overlaying two hardcore dimer configurations gener- 
ates what is known as their transition graph; this graph 
consists of disjoint subgraphs of dimers alternating be- 
tween the two configurations - such subgraphs are shown 
in Fig. g. 

An open graph, as shown on the left of Fig. |l|, must 
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FIG. 1: Transition graphs of initial (black) and final (gray) 
configurations of the pocket algorithm presented in this pa- 
per. Left: In the presence of monomers the graph can be 
open; right: transition loop touching the symmetry axis. The 
cluster move consists in flipping the cluster - Note that we 
may also flip the dimers enclosed by a loop, and that the 
transition graph crosses the symmetry axis at most twice. 



terminate on a monomer. Otherwise, it is a closed loop 
(a "transition graph loop" ) , as the one shown on the right 
of Fig. 0. 

Any Monte Carlo move corresponds to a transition 
graph of the initial and final configurations. The simplest 
transition graphs arise from two pairs of dimers (for the 
square and triangular lattice) or from three (hexagonal 
lattice), which wind around a unit cell of the lattice. A 
Monte Carlo algorithm based on these local moves can 
have the same convergence problems as the local simu- 
lation methods for, e.g., the ferromagnetic Ising model 
near the Curie temperature due to critical slowing down. 
This problem is particularly acute for the square lattice 
dimer model, which has critical correlations. Elementary 
moves with many participating degrees of freedom can 
give rise to a substantial reduction of the Monte Carlo 
correlation time. 

Several algorithms have beep proposed to locate long 
connected transition graphs .Ej Typical problems in this 
context include diverging overhead (most of the time 
spent looking for a long loop is essentially wasted) and 
diminishing return (e.g. a long loop may end up invading 
the full system, or dynamically important moves may not 
be generated for sufficiently large systems). The 'pocket 
algorithm' we present here for hardcore dimers generates 
subgraphs from two configurations related by a global 
lattice symmetry. The algorithm obtains long transition 
graphs with minimal overhead: on lattices of size L x L, 
no moves need to be rejected, and any dimer encoun- 
tered during the construction participates in the move. 
In addition, there is no need for bookkeeping. 

The pocket algorithm proceeds in a series of moves, 
at the beginning of which we randomly pick a symmetry 
axis, and a seed dimer (see Fig. ||). At this initial stage 
of the move construction, the seed is the sole element 
of a set V (the 'pocket', filled with dark dimers), while 
all other (light gray) dimers belong to a set A. Further 
in the move, dimers are shuffled around between sets A 
and V in the following way: At each step, an arbitrary 




FIG. 2: Dimer algorithm, a: Symmetry axis and seed dimer 
(dark gray), sole initial element of the set V, the 'pocket', b: 
An element of the pocket has been reflected with respect to 
the symmetry axis (and transferred from V to A, the set of 
regular dimers). It overlaps with two (dark) dimers, which 
are transferred from A to V ■ The move construction comes 
to an end when V is empty. 

element i of V is reflected with respect to the symmetry 
axis, and moved to A. If i overlaps with other elements of 
A (dark dimers shown in Fig. ||b), the latter are moved 
from A to V. The move is completed as soon as the 
pocket V is empty. During the whole construction, the 
algorithm keeps no memory of the graph, and final suc- 
cess is assured by the underlying symmetry of the lattice 
in combination with the hardcore property of the dimers. 

The pocket algorithm is equivalent to the pivot-cluster 
algorithm,E_3 which has been.-applied successfully to sev- 
eral hard-sphere problems .E30 It provides a greatly sim- 
plified implementation, as the Monte Carlo move is con- 
structed without ever explicitly considering the cluster. 
In the specific application of dimer models, we also use 
reflections with respect to symmetry axes rather than 
point reflections. On the square lattice, for example, 
the algorithm is ergodic if we allow both diagonal and 
horizontal- vertical axes. The first choice allows to change 
the numbers of horizontal/ vertical dimers, the second one 
permits to move through the different winding number 
sectors. The algorithm was tested against exact enumer- 
ation for systems up to size L = 6. 

Transition graphs generated by the pocket algorithm, 
as shown in Fig. [l], are symmetric with respect to the 
symmetry axis, and cross it at most twice. 

An analogous property holds in higher dimensions. 
For dimers on a three-dimensional lattice, for example, 
the appropriate operations are reflections with respect to 
randomly chosen symmetry planes of the (finite) lattice. 
Again, any transition graph is symmetric with respect to 
this symmetry plane, and crosses it at most twice. This 
places a strong constraint on the transition graph, which 
cannot invade space. If we used reflection with respect 
to a point, as was appropriate in other applications ,E£I 
the move would conserve the total number of dimers of 
a given orientation, and also generate frequent transition 
graphs completely invading the lattice. 

Transition loops generated by the algorithm allow to 
flip not only the loop itself but also all the dimers inside 
the (gray) loop area, which is symmetric with respect 
to the axis (see the example in Fig. [I]). This modifica- 
tion, as well as simple generalizations where the pocket is 
initially populated by more than one dimer, or by a spe- 
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cially positioned dimer, have not yet been studied. We 
also note that this algorithm applies to other models that 
can be mapped, onto (not necessarily nearest-neighbor) 
dimer models.E3 

For the undoped system on the square lattice, the 
pocket algorithm generates transition loops of average 
length proportional to L. A more detailed analysis is 
provided in the appendix. 
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III. DIMER MODELS AT FINITE MONOMER 
DENSITY 

We now turn to the problem of the correlations of 
monomers in a background of dimers on the square and 
triangular lattices. We evaluate the monomer two-point 
correlation function for a range of dopings. 

The monomer two-point function is defined as follows. 
Let Z(r) be the partition function of the monomer-dimcr 
problem with a monomer fixed at the origin and another 
at site r: Z(r) counts the number of allowed configura- 
tions subject to the pair of monomers being fixed. The 
two-point function is then defined to be proportional to 
Z(r): C(r) = Z(r)/Z n ; the constant of proportional- 
ity, Zq, is somewhat arbitrary. We choose to normal- 
ize the correlations such that a random distribution of 
monomers would give C(r) — 1. For the two- monomer 
problem however, to make contact with the previous 
literatureJH we normalize the pair correlations such that 
(a) C(l) = l/z and (b) list zC(r) in Tab. @, where z is 
the coordination of the lattice. 

For the square lattice, the two- monomer... correla- 
tions were obtained by Fisher and StephensonQ and by 
HartwigQ by perturbing the Pfafnan matrices introduced 
by Kasteleyn.B For the triangular lattice, there is as 
yet no analytic expression for the asymptotics of the 
monomer correlations J13 and we have obtained the ex- 
act L = 6 results by explicit enumeration. 



TABLE I: Monte Carlo versus exact results for monomer pair 
correlations, zC'(r), on the square (left) and triangular (right) 
lattices. The normalization here is zC(l) = 1. Even distances 
in the square lattice refer to sites on the adjacent row.Q For 
the simulations presented here, the statistical error is about 
10 for both lattices. 
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FIG. 3: Monomer correlations along a coordinate axis: (1) 
pair correlations on the square lattice with L — 64. The fit 
is to an analytic asymptotic expression generalized to include 
periodic boundary conditions. Same-sublattice correlations 
vanish. (2) Four monomers on a square lattice with L — 128. 
The lower curve corresponds to monomers on the same sub- 
lattice, the upper to monomers on different sublattices. The 
upper curve is again compared to the appropriate asymptotic 
function, whereas the lower curve is fitted to a power-law 
r 1,/2 . The deviation of the upper curve from the asymptotic 
expression for large r is due to screening. 



A. The square lattice 

To test the algorithm, first consider two-monomer cor- 
relations, calculated by several authors for the square lat- 
tice. In Fig. H and Tab. [j| we compare our data against 
the analytical results. Note that the fit against the ana- 
lytical asymptotic expression in Fig. || was done using a 
distance variable f = sin (irr / L) / (ir / L) to take into ac- 
count periodic boundary conditions for a torus of length 
LJB 

The fit is against the analytic asymptotic expression, 
and in Tab. ||, we_show a comparison of analytical and 
numerical results E3 The results agree to within the nu- 
merical accuracy of 10 -3 , which improves even further 
for larger monomer densities, p, as the signal-to-noise ra- 
tio increases with the number of monomer pairs as p 2 . A 
more accurate estimate could be obtained by finite-size 
scaling. 



For two monomers on opposite sublattices, the form 
of the correlations, C(r) = a^r^ 1 / 2 defines an en- 
tropic Coulomb interaction between the two monomers: 
V(r) — — ilogr — loga_. Note that whereas the pair 
of monomers are critically confined - C(r) — > alge- 
braically as r — > oo - the expectation value of their sep- 
aration grows with system size: (r) oc L. 

As the square lattice is bipartite, dimers have one end 
on each sublattice. This means that on a torus, there are 
no dimer configurations with the two monomers present 
on the same sublattice: in this case, C(r) = 0. 

The question whether there is an effective interac- 
tion of the same form is nonetheless interesting. To 
avoid the sublattice constraint, we have investigated jour 
monomers on a large lattice of size L — 128. For 
monomers on opposite sublattices, the results are essen- 
tially unaltered from the two monomer case, at least 
for separations which are not too large (r <C L); for 
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monomers on the same sublattice, we now obtain C(r) = 
a + r +1 / 2 , so that the monomers on each sublattice act as 
if they had equal and opposite Coulomb charge. 

The L-dependence of the ratio a_/a + is fixed by the 
fact that the square lattice is at a critical point. This 
allows to write the finite-size scaling ansatz C(r) — 
aL~ 2 c(r J L), where c(x) is a function only depending on 
the ratio of r and L, and a does not depend on L. The 
prefactor L~ 2 follows from the requirement that the area 
underneath the curves remain fixed as L — > oo. With 
c±(x) cx x ±x / 2 , we obtain a_/a + cx L. A fit for L = 128 
and four holes gives a_/a + L/2. 

The dimers thus generate a Coulomb interaction for 
the monomers, V(r) — (qq' /2) logr, with a sublattice- 
dependent sign of the monomer charge q (\q\ = 1). It is 
natural to expect this Coulomb interaction to be screened 
if the monomer density is finite, and their correlations to 
be no longer critical. In fact, for a finite density, p, of 
monomers, a new length scale arises, namely the average 
interparticle spacing, d ~ p^ 1 ^ 2 . 

To analyze the details of the screening analytically, we 
model the monomer system by a classical, continuous 
two-dimensional two-component plasma. With a dimen- 
sionless coupling of (3 = 1/2 (the prefactor of the loga- 
rithm) , this plasma is located well in its high-temperature 
phase. This phase extends up to (3 = 2, where the contin- 
uum approximation breaks down as the plasma collapses 
in the absence of a repulsive hard core to the interaction. 

In this high-temperature regime, we expect Debye 
screening to be operative. The resulting long-distance 
connected correlations have been obtained in Ref. |22| 
(where an expression going beyond Debye screening is 
given); they are 



C(r) = 1 ± X 
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with the Debye screening length: Ad = l/^2nf3p, and 
with x = P to lowest order in the coupling. The choice 
of sign depends on the sublattice. 

With the new length scale, Ad cx d, the criticality of 
the correlations is destroyed, although the critical be- 
havior remains visible for r <C Ad- This is why it was 
possible to find the two-particle interaction from the four- 
monomer problem above. 

In order to determine whether this scenario is indeed 
realized for the monomer-dimer mixture, we have plot- 
ted C(r) versus the scaled distance r/Xo for varying 
monomer densities for a system of size L = 64 (Fig. [|). 

From the form of the expression for C(r), all curves 
should collapse on top of one another. However, for a 
finite system, one has to take into account the presence 
of two length scales, L and Ad, with the added feature 
that Ad < L, so that a simple one-parameter finite-size 
scaling ansatz no longer applies. For Ad < L (large 
monomer densities) , the collapse works indeed very well. 

The data do not collapse at small densities, which 
is a finite-size effect: In Fig. |^, we plot the monomer 
correlations at a fixed monomer density p = 2/64 2 for 
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FIG. 4: Monomer correlations C for different monomer densi- 
ties, against the scaled distance variable r/Au, for the square 
lattice with L — 64. 
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FIG. 5: Monomer correlations C against the scaled distance 
variable r/Ao at fixed monomer density p = 2/64 2 for square 
lattices of size L = 64, 128, 256. 



L = 64, 128 and 256. Whereas the two-hole problem 
suffers from the sublattice pathology, one can see that 
the data for L = 128 and 256 collapse nicely for short 
distances, with disagreement arising for distances around 
L 1 2 where the periodic boundary conditions become im- 
portant. 

In the region where they collapse, the curves in Figs. ^ 
and U agree with the asymptotic Debye prediction for 
large values of r/\r>. At small r, the effectively un- 
screened behavior is visible, leading again to a unique 
curve for C(r). 



B. The triangular lattice 



Dimer correlatior 



to be short-rangec 
the square lattice! 



_on the triangular lattice are known 
II and not critical as is the case for 
This leads one to expect that the 
monomer correlations will also be short-ranged as the 
presence of a monomer should have little effect at a dis- 
tance large compared to the dimer correlation length. 

We have plotted the correlations of a pair of monomers 
in the inset to Fig. |^. The correlation length indeed is so 
small that already for a triangular lattice of size L = 12, 
finite-size effects are no longer detectable. This is also 
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FIG. 6: Monomer correlations C, along a coordinate axis at 
variable doping for a lattice of size L = 24. The abscissa is in 
units of lattice spacings. Correlations are normalized against 
a random distribution. Inset: two-monomer correlation in 
a lattice of size L — 12, In both cases, the a>axis denotes 
separation along one of the lattice vectors. Note that the 
y-axis only varies by less than 15%. 



shown in Tab. |. 

At 'large' distances, the correlations tend to a constant, 
nonzero value of C(oo) = 0.1495 ± 0.0001. (This value is 
given, in keeping with the convention of the square lat- 
tice case, with a normalization such that C(l) = 1/6.) 
They do so with an extremely short correlation length, 
which we cannot determine accurately as there are only 
about 3 data points above the noise level. This makes 
it hard to fit an asymptotic form, as the presence of a 
power-law and/or oscillatory multiplicative factor can- 
not be detected on this basis. From the absence of any 
signal for distances beyond r = 3, we conclude that the 
correlation length is less than one lattice spacing. 

This nonzero large-distance limit of the monomer cor- 
relation function implies that monomers in the classical 
hardcore dimer model on the triangular lattice are decon- 
fined - the free energy for monomers a distance r apart 
does not diverge as r — > oo. This contrasts with the 
square lattice, where this quantity diverges logarithmi- 
cally. 

Upon further doping, the correlations barely change 
(see Fig. ^) ; a small decrease of the correlations only be- 
comes visible for a doping level of 1/8. This is not unex- 
pected as the average monomer spacing always remains 
well above the two-monomer correlation length. 



IV. CONCLUSION 

We have studied the correlations of monomer-dimer 
mixtures on the square and triangular lattices with a 
new Monte Carlo algorithm. The monomers interact via 
the background of dimers: the critical correlations of the 
dimers on the square lattice give rise to a two-monomer 
interaction of long-range Coulomb nature, whereas the 
disordered triangular correlations generate an exponen- 
tially decaying one. 

The triangular model is deep in its disordered phase, 



FIG. 7: Distribution of transition loop lengths for the square 
lattice with L — 64. 
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FIG. 8: Distribution of transition loop lengths for the trian- 
gular lattice with L — 64. 



where it remains in the presence of a finite density of 
monomers. The square model, located at a phase bound- 
ary, is driven into a disordered phase on doping. This ef- 
fect is rendered analytically by modeling the monomers 
as a two-dimensional two-component plasma, with the 
long-range correlations given by a simple Debye-screened 
form. 



Appendix 

In this appendix, we analyze the distribution of tran- 
sition graph lengths in the undoped dimer model, where 
only loops are generated. These loops can be of several 
types: disconnected pairs of loops not touching the sym- 
metry axis, connected loops passing through the symme- 
try axis, but not winding around the box, and finally 
those winding around the box. Each of these distribu- 
tions strongly depends on the lattice and on the details 
of the choice of axes. The distribution of loop lengths, A, 
generated by the pocket algorithm is shown in Fig. ^ for 
the square lattice. Whereas the disconnected loops fol- 
low a nice power law distribution P (A) disconnect ~ A 5 
for A < L, and are always made up of an even num- 
ber of dimers, the connected loops are distributed as 
P(A) conne ct ~ A -075 . For the reflection about a di- 
agonal, the loop lengths are always even, whereas they 
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are odd for the horizontal- vertical reflections that do not 
wind around the lattice. If a loop does wind around the 
box, its length can be odd or even. Notwithstanding this 
intricate structure, we find that the mean loop length 
grows as L^, with fj, = 1.00 ± 0.01. On the triangu- 
lar lattice, the distribution of loop lengths is also nicely 
split into connected and disconnected loops (see Fig. ||) . 
There, we find for the mean loop length /i = 1.47 ± 0.02. 
We note that the algorithm is biased towards longer sub- 
graphs of the transition graph of a configuration with its 
reflection with a factor proportional to their length. 
Finally, it would be interesting to analyze the expo- 



nents in more detail, in particular as the square model 
maps onto-,a height modelJ£3 whereas the triangular one 
does not .Ell 
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